function diff = transition_t_OJS(theta, J_old, Q_ini, param, aux, ind_min)

    G_r         = nan(size(aux.t_vec));
    hires_tot   = nan(size(aux.t_vec));
    new_chains  = nan(size(aux.t_vec));
    mean        = nan(size(aux.t_vec));
    mean_m      = nan(size(aux.t_vec));
    mean_w      = nan(size(aux.t_vec));
    mean_lnw    = nan(size(aux.t_vec));

    excess=nan(size(aux.t_vec));
    m_l=nan(size(aux.t_vec));
    m_h=nan(size(aux.t_vec));
    m_u=nan(size(aux.t_vec));

    u_t   = nan(size(aux.t_vec));
    ete   = nan(size(aux.t_vec));

    %% HJB
    for t=length(aux.t_vec):-1:1

        if t==1
            aux.dt=aux.t_vec(t);
            param.mu = param.mu_base+(aux.p_t(1)-(1-param.shock))/aux.dt;
        else
            aux.dt=aux.t_vec(t)-aux.t_vec(t-1);
            param.mu=param.mu_base+(aux.p_t(t)-aux.p_t(t-1))/aux.dt;
        end

        param.lambda=max(0,interp1(aux.t_lam,theta,aux.t_vec(t),'linear','extrap'));

        J_i=HJB_t_OJS(J_old,param,aux);

        J_old=J_i;

        m_l(t)=interp1(J_old,exp(param.ln_m_HJB),0,'pchip');
        m_h(t)=interp1(J_old,exp(param.ln_m_HJB),param.c,'pchip');

        J_old=min(param.c,max(0,J_old));

     end

    %% FPE

    for t=1:1:length(aux.t_vec)

        % Set time specific variables 
        param.lambda=max(0,interp1(aux.t_lam,theta,aux.t_vec(t),'linear','extrap'));
        param.m_l=interp1( aux.t_vec, m_l,aux.t_vec(t),'pchip');
        param.m_h=interp1( aux.t_vec, m_h,aux.t_vec(t),'pchip');    

        x0= fzero(@ (x) delta_OJS(abs(x)+param.m_h,param), [0,1], optimset('display','off'));
        param.m_u=param.m_h+abs(x0);
        m_u(t)=param.m_u;     


        if t == 1

            Q_old= Q_ini_fun_OJS(Q_ini,param,aux);
            m_old=[param.m_l, param.m_u];

            ln_m_ini = ln_m_KFE(m_old,1,aux);
            q_ini=(Q_old(2:end)-Q_old(1:end-1))./(ln_m_ini(2:end)-ln_m_ini(1:end-1));

            mean_ini =((q_ini.*(ln_m_ini(2:end)-ln_m_ini(1:end-1)))'*exp(ln_m_ini(2:end)./(1-param.alpha)));

            aux.n=50;
            aux.dt=aux.t_vec(t)/aux.n;
            param.mu = param.mu_base+(aux.p_t(1)-(1-param.shock))/(aux.n*aux.dt);
            Q_old=KFE_t(m_old,Q_old,delta_fun(exp(ln_m_KFE([param.m_l,param.m_u],1,aux)),param),param,aux);        
            aux.dt=aux.t_vec(t).*aux.n;

        else        
            aux.n=1;
            aux.dt=(aux.t_vec(t)-aux.t_vec(t-1))./aux.n;
            param.mu=param.mu_base+(aux.p_t(t)-aux.p_t(t-1))/aux.dt;

            Q_old=KFE_t(m_old,Q_old,delta_fun(exp(ln_m_KFE([param.m_l,param.m_u],1,aux)),param),param,aux);
        end 

        m_old=[param.m_l, param.m_u];

        ln_m=ln_m_KFE([param.m_l,param.m_u],1,aux);

        q=(Q_old(2:end)-Q_old(1:end-1))./(ln_m(2:end)-ln_m(1:end-1));
        mean_m(t) = mean_t_fun(exp(ln_m), Q_old);
        mean_w(t) = mean_t_fun(wage_fun(exp(ln_m),param), Q_old);
        mean_lnw(t) = mean_t_fun(log(wage_fun(exp(ln_m),param)), Q_old);
        ete(t)  = jtj_num(delta_fun(exp(ln_m),param), Q_old);
        u_t(t)    = Q_old(1).*param.s;
        G_r(t)          = interp1( exp(ln_m), (Q_old - Q_old(1))./(Q_old(end) - Q_old(1)), param.m_h, 'linear');
        hires_tot(t)    = ete(t).*(1-u_t(t))+param.lambda*u_t(t);
        new_chains(t)   = param.lambda*u_t(t) + param.s.*param.lambda.*G_r(t).*(1-u_t(t));

        mean(t)=((q.*(ln_m(2:end)-ln_m(1:end-1)))'*exp(ln_m(2:end)./(1-param.alpha)));

        if t==1
            excess(1)=-(mean(t)-mean_ini) - (mean(t)-mean_ini)/20;
        else
            excess(t)=-(mean(t)-mean(t-1)) - (mean(t)-mean_ini)/20;
        end
        excess(t)=excess(t)*30/mean_ini;

        if ind_min == 0 
            if t == length(aux.t_vec)
                diff.F_p_t( end+1, :) = -delta_p_fun(param.m_grid,param)./param.s./param.lambda;

                diff.F_t( end+1, :) = 1-delta_fun(param.m_grid,param)./param.s./param.lambda;
                diff.G_t( end+1, :) = interp1( [0; exp(ln_m(1))/2; exp(ln_m); exp(ln_m(end))+1;exp(ln_m(end))+10], [0; 0; (Q_old-Q_old(1))./(Q_old(end)-Q_old(1)); 1; 1], param.m_grid);
                diff.t_dist( end+1) = aux.t_vec(t);
            elseif (aux.t_vec(t) < t_grid_dist(ind_del) && aux.t_vec(t+1) > t_grid_dist(ind_del))              
                diff.F_p_t( end+1, :) = -delta_p_fun(param.m_grid,param)./param.s./param.lambda;

                diff.F_t( end+1, :) = 1-delta_fun(param.m_grid,param)./param.s./param.lambda;
                diff.G_t( end+1, :) = interp1( [0; exp(ln_m(1))/2; exp(ln_m); exp(ln_m(end))+1;exp(ln_m(end))+10], [0; 0; (Q_old-Q_old(1))./(Q_old(end)-Q_old(1)); 1; 1], param.m_grid);
                diff.t_dist( end+1) = aux.t_vec(t);            
            elseif (aux.t_vec(t) > t_grid_dist(ind_del) && aux.t_vec(t-1) < t_grid_dist(ind_del)) || aux.t_vec(t) == t_grid_dist(ind_del)             
                diff.F_p_t( end+1, :) = -delta_p_fun(param.m_grid,param)./param.s./param.lambda;

                diff.F_t( end+1, :) = 1-delta_fun(param.m_grid,param)./param.s./param.lambda;
                diff.G_t( end+1, :) = interp1( [0; exp(ln_m(1))/2; exp(ln_m); exp(ln_m(end))+1;exp(ln_m(end))+10], [0; 0; (Q_old-Q_old(1))./(Q_old(end)-Q_old(1)); 1; 1], param.m_grid);
                diff.t_dist( end+1) = aux.t_vec(t);            
                ind_del = ind_del + 1;
            end
        end
    end

    if ind_min == 1
        diff.excess =  interp1(aux.t_vec,excess,aux.t_lam, 'pchip','extrap');
        diff.mean=mean;

        diff.mean_m=mean_m;
        diff.ete = ete;
        diff.u = u_t;

        diff.lambda=theta;
        diff.t_grid= aux.t_vec;

        diff.hires_tot= hires_tot;
        diff.new_chains= new_chains;

        diff.m_l_t=m_l;
        diff.m_h_t=m_h;
        diff.m_u_t=m_u;

    elseif ind_min == 0 
        diff.u      = interp1(aux.t_vec,u,diff.t_grid,'pchip','extrap');
        diff.m_l    = interp1(t_grid_l,m_l_t,diff.t_grid,'pchip','extrap');
        diff.m_h    = interp1(t_grid_m,m_h_t,diff.t_grid,'pchip','extrap');
        diff.m_u    = interp1(aux.t_vec, m_u, diff.t_grid,'pchip','extrap');
        diff.lambda = interp1(aux.t_vec, theta, diff.t_grid, 'pchip', 'extrap');
        diff.zeta   = param.zeta.*ones(1, N);
        diff.ete    = interp1(aux.t_vec, ete, diff.t_grid, 'pchip', 'extrap');
        diff.mean   = interp1(aux.t_vec, mean, diff.t_grid, 'pchip', 'extrap');
        diff.mean_m = interp1(aux.t_vec, mean_m, diff.t_grid, 'pchip', 'extrap');
        diff.mean_lnw = interp1(aux.t_vec, mean_lnw, diff.t_grid, 'pchip', 'extrap');
        diff.mean_w = interp1(aux.t_vec, mean_w, diff.t_grid, 'pchip', 'extrap');
        diff.p      = interp1(aux.t_vec, aux.p_t, diff.t_grid, 'pchip', 'extrap');
        diff.hires_tot= hires_tot;
        diff.new_chains= new_chains;
    end
    
end